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Abstract: We report calculations of the parameter Bk appearing in the AS" = 2 
neutral kaon mixing matrix element, whose uncertainty limits the power of unitarity 
triangle constraints for testing the standard model or looking for new physics. We 
use two flavours of dynamical clover- improved Wilson lattice fermions and look for 
dependence on the dynamical quark mass at fixed lattice spacing. We see some 
evidence for dynamical quark effects and in particular Bk decreases as the sea quark 
masses are reduced towards the up/down quark mass. 
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1. Introduction 

Bk is the AS = 2 neutral kaon mixing matrix element normalised by its vacuum 
saturation approximation (VSA) value, 

i%(^) - 8 f2yrl 2 ' ( L1 ) 

3 JK' n K 

with /i indicating the scale dependence of the operator Q(p) = S7 M (1 — 75)^ S7 M (1 — 
75)0?. This can be related to the one-loop renormalisation group invariant (RGI) 
value Bk through 
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where n/ is the number of active flavours at the relevant scale, and 70 and (3q have 
the scheme independent values of 4 and 11 — 2n//3. We use the MS scheme for which 
J is calculated to NLO in [1]. To go from MS at 2 GeV to the RGI value we note 

that there are four active flavours and ^ms(^) = 1-792. Starting from the PDG value 

(5) 
of Aq^, d = 216 MeV and matching the strong coupling at the charm threshold we 

obtain Bk = 1.404 Bk(MS, 2 GeV). This is rather insensitive to the value of n/ [2]. 

The standard model expression for the indirect CP violating parameter [3] as 

quoted in [4], 

e K = rjA 2 B K [1.11(5) • A 2 {1 - p) + 0.31(5)] , (1.3) 

defines a hyperbola in the (p, fj) plane, A, p and fj being parameters of the CKM 
matrix elements and unitarity triangle [5]. The theoretical uncertainty in the value 
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Table 1: Some previous lattice calculations of Bk- NP refers to non-perturbative renor- 
malisation. Only the last number is unquenched. 

of Bk remains the dominant uncertainty when we try to use this expression along 
with the experimental value of Ek to constrain the triangle. This has resulted in a 
great deal of activity in the lattice community to refine this calculation. 

There is a relatively long history of Bk calculations in different frameworks. 
Some recent lattice calculations are listed in table 1. A more comprehensive sum- 
mary can be found in [16], with numbers from other methods dispersed over a rela- 
tively wide range. Over the years the quenched lattice value of Bk has more or less 
settled down. The 1997 quenched value of .Bj^MS, 2 GeV) = 0.63(4), corresponding 
to Bk = 0.87(6), using staggered fermions [7] remains the benchmark and is the 
value usually quoted for phenomenology. Other quenched numbers are more or less 
consistent with this. The error quoted however does not include any estimate for 
quenching effects and this has been estimated to be as high as 15% [17]. Unquenching 
remains the primary systematic effect to be addressed. 

There has been one preliminary report of a complete unquenched calculation us- 
ing Domain Wall (DW) fermions from the RBC collaboration [15] and a few other at- 
tempts on selected sets of parameters using Wilson and staggered fermions. Though 
the central values for Bk from DW fermions have often been on the lower side, the 
unquenched DW preliminary number is really at the lower end of the spectrum. 
Other attempts to unquench, e.g. [4, 18, 19, 20, 21, 22], have not always been able 
to see a definite effect. However, it has been noted [23] that though the unquenched 
numbers are consistent with the quenched numbers within errors, they are system- 
atically lower. Hence, it is difficult to reach an unambiguous conclusion on the true 
effect of dynamical fermions on this quantity. 

In this paper, we report on a calculation using two degenerate flavours of dynam- 
ical (clover-improved) Wilson fermions. In order to look for sea-quark dependence in 
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Bk we use three different sea quark masses in the region rap j my > 0.7 on a volume 
of 16 3 x 32 (mpL > 7) but a nearly constant lattice spacing. To achieve the latter a 
set of values of the bare coupling and bare dynamical quark mass have been chosen 
in [24, 25] to keep the lattice spacings, defined using the scale r [26], as fixed as 
possible. 

For Bk, we see some evidence for dynamical quark effects and in particular the 
values decrease as the sea quark mass decreases from the simulated range towards 
the up/down quark values. 

This calculation is undertaken as an intermediate step towards a complete un- 
quenched evaluation of Br. In the near future one might hope to perform detailed 
studies over lighter and larger samples of sea quark masses at different lattice spac- 
ings in order to make the continuum extrapolation. In the meantime, exploratory 
studies may help as a guide to those regions of parameters accessible today. 

The plan of this paper is as follows. In section 2 we give the basic definitions 
and in section 3 introduce the quantities relevant for a lattice estimate of Bk- In 
section 4, we discuss the analysis and present our values and then we have some 
concluding remarks in the final section. 

2. Setup of the calculation 

In the continuum, the operator of interest in eq. (1.1) is 

Q AS=2 (fi) = QM = s^d s-fd + s^ l5 d s^lsd, (2.1) 

which is the parity conserving part of Q(n) in eq. (1.1). For Wilson fermions, owing 
to the breaking of explicit chiral symmetry, there is a mixing of this operator with 
other four-fermion operators. Therefore one has to work with a complete basis of 
operators and subtract the extra ones. One such set is 

Qi(li) = s'j^d s'fd + s'j^d s7 M 7 5 d 

Q 2 (aO = si»d sj^d - 57^75^ s^j 5 d 

Qz(fi) = sdsd + S'-y^d s^d (2.2) 

Qa{ij) = sdsd — S75G? S75G? 

Qb(lj) = sa^dsa^d. 

Together with the overall multiplicative renormalisation, the subtraction of the un- 
wanted operators may be expressed in a compact form as 

Q nM (t<) = Z(,,.<fi) [Qr + yiMyfiCft*]. (2.3) 
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The renormalisation coefficients Z and Aj have been determined perturbatively 
for MS-NDR in [27, 28]. Once the renormalisation and subtraction of eq. (2.3) is 
carried through, we have the matrix element for our desired operator in eq. (1.1). 

For fermion implementations which (nearly) respect chiral symmetry, e.g. in [6, 
11, 13], the chiral behaviour is not modified by lattice artefacts and Bk{p) can be 
obtained from matrix elements of kaons at rest. But for Wilson fermions as, for ex- 
ample, in [4, 18], lattice artefacts introduce chiral symmetry breaking contributions 
to Bk in the chiral limit. In our case, even though we use an improved-clover action, 
four- fermion operators are unimproved and 0(a) artefacts may be present. To par- 
tially remove them at finite lattice spacing another degree of freedom is required and 
this can be done by introducing non-zero momentum kaons. Simulations at different 
lattice spacings and extrapolation to the continuum also allow lattice artefacts to be 
removed. 

Let us now consider matrix elements with non-vanishing external momenta and 
generic pseudoscalar mesons. On the lattice, the chiral behaviour of the matrix 
element with non- vanishing external momenta can be parametrised as [17] 

(P°, p |Q(/i) \P°, q) =a' + (5'rnp + 8'mp + 

(P ■ q) (7 + i + (e + e')m 2 p + (f + 0(p •?)) + ••• (2-4) 

where all the quantities are expressed in lattice units and the ellipsis stands for higher- 
order terms in p-q and mp. All the primed coefficients are lattice artefacts. However, 
while 7' and e' are corrections of (9(a) to the corresponding physical contributions, the 
parameters a', 13' and 8' are absent in the continuum limit and have to be subtracted 
from the estimate of Bk in eq. (1.1). In particular the a' term makes Bk divergent 
in the chiral limit. 

For our calculation with Wilson fermions, we neglect higher order terms and use 
the following expression for the matrix elements: 

(P ,p\Q(fi)\P ,q)=a' + [3'm 2 P +( 1 + 1 ')(p-q). (2.5) 

3. Numerical simulation 

In this work Bk is calculated using Clover-improved Wilson fermions [29] on the 
UKQCD set of unquenched configurations listed in table 2. Details of the generation 
of the gauge configurations can be found in [24, 25]. To have a decorrelated sample, 
configurations separated by 40/50 trajectory steps are used. These configurations 
are on a lattice of 32 x 16 3 points. 

The lattice spacings determined from the Sommer scale, ro, are very similar 
for these sets. However, there are concerns that the /t sea -dependence of the lattice 
spacing observed in these configurations is due to the proximity to a phase transition 
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Set (3 c S w ^sca a(fm)[GeV~ ] (mp/m v ) KB ^ =Kval No. of configs 

I 5.20 2.0171 0.1350 0.103(2) [1.91(2)] 0.70(1) 100 

II 5.26 1.9497 0.1345 0.104(1) [1.90(2)] 0.78(1) 100 

III 5.29 1.9192 0.1340 0.102(2) [1.94(2)] 0.83(1) 80 

Table 2: The configurations used. Values for lattice spacings are as calculated from the 
value of the scale, ro, in lattice units from the UKQCD set [24, 25]. 

around a ~ 0.1 fm where there may be large cutoff effects in the dynamical case 
[30, 31, 32] and therefore needs to considered with caution. We take the view that, 
nevertheless, these sets of configurations do have some degree of matching according 
to a valence-quark-independent definition of an effective lattice spacing, and thus, 
unless our physics is completely overwhelmed by any nearby phase transition, a 
combined analysis of the data as a function of /t sca is worthwhile. It may be noted 
that since Bk is dimensionless, the lattice spacing enters through discretisation errors 
but not via an overall power of a. Moreover, when analysing the sea quark mass 
dependence, we use the variable (amp) 2 (ft sea , rc sea ) which in our case is equivalent to 
using (r mp) 2 since our lattice spacing is defined through r and r^vrip = (r /a)xamp 
with (r /a) fixed for these lattices [24, 25]. 

Propagators and correlators were calculated using the FermiQCD [33, 34] code. 
Five valence quark propagators at k — 0.1356, 0.1350, 0.1345, 0.1340 and 0.1335 
were generated for each sea quark using the Stabilised Biconjugate Gradient method 
[35]. Smearing was tried, but since it did not give any significant improvement in 
the signal, the results presented here are for the local case (see comments in the next 
section). 

To calculate the matrix element on the lattice the standard procedure [36] is 
followed where we calculate the 3- and 2-point correlation functions 

C®\t x ,ty,p x , P y,n) = ^(^(^^(O^jPsi^gje^e^ (3.1) 

xy 
-t yJx >>0 Zpe -E x t x ^ Q ^ p)Zpe -E v t y ^ 

CfjfrP*) = £U(*,*).#,0)) e ** (3.2) 



**>P -2- -2- f ~E x t 



ZjiZj.e 



Here the J's are kaon interpolating operators and can be of pseudoscalar or axial 
vector type; they can also be local or smeared. We use pseudoscalar current sources. 
In the 3-pt functions, the operator is fixed at the origin and t y is kept fixed at a 
particular value, while t x is varied over the full temporal range of the lattice. The 
main reported results are for t y = 10. We have checked with other neighbouring 
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Figure 1: Fits for lattice matrix elements for the complete set of bare operators for a 
sample of our data (set I, K va i = 0.1350). Ratios of the 3-pt correlators to two 2-pt (PP) 
correlators are fitted in the interval t x = 22 — 27 for t y = 10 (see eq. 3.4). Correlators are 
shown for zero momentum. The fitted ones are those of interest (P \Qi\P ) while the other 
plateau in the first half of the lattice corresponds to the (P P \Qi\0) matrix elements. 

values of t y but observe no dependence, implying that the ground state is reason- 
ably well isolated by this time. For the momentum configurations, we have chosen 
{Px,p y } = {(0,0,0), (0,0,0)}, {(0,0,0), (1,0,0)} and {(1, 0, 0), (0, 0, 0)} where the 
average over equivalent configurations is understood. 
For simulation, we use the simpler basis of 

Qv(fj) = sj^dsj^d 

Qa{h) = ~sivl<r,d~sYlbd 

Qi(n) = sd sd (3.3) 

Qp(h) = s^ 5 dsj 5 d 

Qs(fj) = sa^dsa^d, 

which is related to our renormalisation basis introduced in the previous section 
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Table 3: Perturbative matching coefficients to go from B^^fi 
GeV). 



I /a) to B% b (n = 2 



through a simple rotation. Fitted ratios for this basis that give us the matrix el- 
ements in lattice units, Q' att , (i = V, A, I, P, S) are plotted in fig. 1. 

To go directly to MS at ji = 2 GeV, we note that in our case (a/x) ~ 1 and we 
can naively use standard perturbation theory at one-loop. For the coupling there is a 
range of choices that may lead to different numerical values. We use the boosted bare 
lattice coupling, g$ = 6//5(P), where (P) is value of the relevant average plaquette 
and our values are {0.5336, 0.5399, 0.5424}. For csw the one-loop value of 1.0 is 
used. The perturbative matching coefficients thus obtained are listed in table 3. 

To extract the desired matrix element the following ratios are formed: 
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where Z& is the axial current renormalisation. 
At this stage one may fit the equation 
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to obtain estimates for P>k from 7 [8, 37], by neglecting 7'. In the fit, the parameters 
with tildes are taken to be constant and hence the estimates are for effective values 
of Z-p and fp in our range of simulation. In this manner, for a set of different valence 
quarks with a given sea quark mass, this approach gives an estimate of the leading 
term in an expansion of Br for that set with the kaons not necessarily being at the 
physical kaon mass. 

To obtain estimates of Br for each (« sea , «; va i) combination, which will then 
allow us to extrapolate in the quark masses, we follow the approach of [18, 27]. 
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Let us call the non-zero- and zero-momentum R^s Rzip) and i?3(0) respectively, 
corresponding to X(p) and X(0) defined in eq. 3.4. The two non-zero momenta 
{Px,P y } = {(0, 0, 0), (1, 0, 0)} and {(1, 0, 0), (0, 0, 0)}, have been averaged, since they 
are estimates of the same matrix elements in the continuum and indeed numerically 
are found to be very similar. Then we have 



R 3 {p)-Rs(0) 



X(p)-X(0) 



— Bk((1, Ksea, ^val)- (3.£ 



These can then be used in our chiral extrapolations in the sea and valence quarks. 
At the same time, fitting these values to a constant for a given sea quark is similar 
to estimating 7 from a fit of eq. 3.7. At higher orders of momentum, this expression 
differs from the correct dependence of Bk by a term like C,mpE(p) [27]. We have 
found the coefficient £ of this term difficult to determine, particularly for our limited 
set of momenta. However, if we were able to make this correction, it would simply 
change our values of Bk within our systematics, leaving our conclusions unchanged. 

4. Analysis and discussion 

The values obtained for Bk(M.S,2 GeV) for our sets of masses are tabulated in 
table 4. We refer to the ones quoted from eq. (3.7) following [8, 37] and from eq. (3.8) 
following [18, 27] as method I and II respectively. 

We have degenerate valence quarks. So, SU(3) breaking effects due to m s 7^ 
m U( i are not accounted for. Rather our kaon is made up of two quarks around 
m s /2. Moreover, the results in table 4 are obtained for local sources. Indeed, we 
have not seen any significant improvement of the signal from smearing. This is 
not unexpected since we have a local operator at the origin and can smear only at 
the sink, which is usually less effective than source smearing. It may also be due 
to a lack of optimisation of the smearing parameters. However, results were fully 
compatible with those using local operators and we have restricted the presentation 
to the simpler case. 

In fig. 2, we plot -Br-(MS, 2 GeV) from method II as a function of the correspond- 
ing squared pseudoscalar masses over the complete set of our valence and sea quark 
masses. We observe the points for the lightest valence quarks diverging for the differ- 
ent sea quarks. Here, it may be noted that for K va \ ^> n sea the theory becomes more 
like quenched. This effect is clearly seen in fig. 3 of [38] where as the valence quark 
becomes lighter the partially quenched curves leave the full theory and tend towards 
the quenched one. Finite volume effects are, in general, expected to be small [39], 
but for some regions of parameter space, particularly for very light quarks, it has 
been suggested that finite volume effects can obscure the chiral behaviour in Bk [38]. 
The JLQCD collaboration [40] observes finite volume effects for lighter sea quarks 
for the same action, but for our parameters they have excluded finite volume effects 
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Table 4: Simulated values of B K (MS, 2 GeV). Method I refers to a direct fit of eq. (3.7); 
while in method II, eq. (3.8) is used to obtain values for each (n sca , n va \) combination. 



for pseudoscalar meson correlators down to just beyond our lightest point in set I. 
Indeed we find the finite volume correction from [38] to be —0.1% for this point. 
Nonetheless, we note that, contrary to the other sets, for set I, the 0(a) discretisa- 
tion error parameters a and (3 turn out to have finite values of —0.06(2) and 0.23(8). 
The effects of these terms are greater at lighter quark masses and we cannot be sure 
that the curvature observed here is due to a true chiral behaviour. As can be seen 
from our values of rap/my, this point is at a considerably lighter mass than all the 
other points. Therefore, we choose to be cautious and exclude it from our analysis. 
It would be interesting to know if non-perturbative renormalisation [41, 42], and/or 
the improvement programme of [43, 44] could lead to better chiral behaviour. 

Now let us consider the values from method I. It is notable that for these rather 
heavy sea quarks these numbers are compatible with quenched estimates. This is the 
reason that previous attempts to unquench for a fixed heavy sea quark mass have 
found it difficult to disentangle the unquenching effects. 

Since we have more than one sea quark mass in our simulation, we can attempt 
to extrapolate these numbers to the chiral limit. We use a linear fit versus the unitary 



pseudoscalar masses (amp) 



IK K 



K, va \) and go to the up/down limit. This gives us 



B K (MS, 2 GeV) = 0.49(13), 
which corresponds to Bk = 0.69(18). 
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Figure 2: Values of Bk(MS, 2 GeV) for each (K sca , n va \) combination plotted as a function 
of the corresponding squared pseudoscalar masses. The dashed lines joining the points are 
just for a visual guide separating the sets with different sea quarks. The filled points joined 
by a solid line are the unitary ones for which K sea = K v3b \. The lightest point for set I 
(marked by a large cross) is excluded from the analysis. 



In this method we estimate 7 in eq. 3.7. As mentioned in the previous section, 
the valence quarks are not necessarily such that mp = m? K yh . In fact one can note 
by comparing with the last column of table 4 that these values are in the simulated 
region. Therefore one may think of this estimate as one of Bk where the sea quarks 
are realistically light but the valence quarks are heavier than the physical strange 
quark. 

A somewhat complementary approach, would be to follow the route of [15] and 
take the unitary points, i.e. the points with K sea = K va \ from method II, for extrapo- 
lation to the physical kaon mass [fig. 4]. This leads to 



B K (MS, 2 GeV) = 0.48(13), 



(4.2) 



Corresponding to Bk = 0.67(18). Here we have a more reasonable valence mp 



m 



phys 
K 



but on the other hand the sea and valence quarks are degenerate and hence 
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Figure 3: Fit to the data from method I. The values quoted is from the linear extrap- 
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the sea content is not as light as the up/down quarks. To understand how much this 
may affect us we note that if we take all the quark masses (both valence and sea) to 
zero our value of Bk goes down to 0.40(17) and Bk = 0.56(24). 

A combined analysis of valence and sea quarks has been tried for the spectroscopy 
studies in [25, 40]. With more momenta, higher statistics and/or a larger sample of 
sea and valence quark masses and if the higher order terms in Bk could be estimated, 
this would be a possible route to an estimate of Bk at the physical valence and sea 
masses. 

Even though we recognise that the presence of several artefacts does not allow a 
quantitative estimate of the sea quark dependence, it does seem that dynamical quark 
effects can be quite significant. There also seem to be indications that, incorporating 
dynamical quarks lowers the value of Bk- Taking this together with the observation 
in [23] that the Nf = 2 numbers are always lower than those for Nf = 0, a statement 
also valid for subsequent works, we see that when one has two finite-mass but still 
heavy sea quarks, Bk starts to decrease but is still consistent with the quenched 
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Method II 
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Figure 4: Unitary fit of the data. The value quoted is from the linear extrapolation, 
whereas the quadratic and chiral log-type fits are added for illustration. The vertical lines 
show the positions where mp = rn^ ys and mp = m K ys . respectively. 



value within errors. When the sea quarks can be taken to the massless limit, the 
value of Bk becomes distinctly lower than the quenched result. It is also intriguing 
to note that in a recent study where Bk is taken as a free parameter and fitted using 
the other unitarity triangle constraints, the value obtained is Bk = 0.69(11) [45], 
again lower than the usual quenched lattice value. 

Owing to the exploratory nature of our analysis and large statistical errors, a 
study of systematic errors such as those connected to choices of fit window, chiral ex- 
trapolation, renormalisation method, the fixed time at one end, the strong coupling, 
Aqcd, etc. has not been addressed. 

5. Conclusion 

We have presented results for Bk calculated using non-perturbatively (9(a)-improved 
Wilson fermions with two dynamical flavours for three sets of relatively small vol- 
ume lattices of matched spacing. Despite some concern about the robustness of the 
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estimates due to various lattice uncertainties, there are indications that dynamical 
quark effects are important and lead to a lower value of Bk- 
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